%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Of Cities and Slums: Main Code for Solving Equilibrium of the Extended 
% Version of the Model, with a Fourth Location: Distant Poor Neighborhood
% % ---- Calls the Functions: 
%     RuralEquilibriumFunction.m,
%     M_OnlyFunction.m and 
%     UrbanEquilibriumFunction.m
%% Outline
% 
% I.      Setting Up the Economy: Parameters and Functional Forms
% II.     Setting Up the Code: Computational Aspects, e.g. gridpoints, and options.
% III.    Initial Conditions: Values Set for population locations, etc
% IV.     Big Loop: Equilibrium location decisions, q_R, q_F, q_C
%             Inner Loop: Equilibrium Prices and Quantities, given locations.
% V.      Reporting the Results      


clear all
close all
warning off


%% Preferences parameters
   %--Stone-Geary for Consumption
        alphaS = 0.147;
        alphaA = 0.01;
        alphaM = 1 - alphaA - alphaS;
        cSbar  = 0.2;
        cAbar  = 1.9; 
   %--Curvature in the utility of Consumption   
        gam = 1;
   %--Idiosyncratic Locational Preferences   
        etaE   =  1.889;
        etaE_R =  1.889;
   %-- Utility Weight of Children's Education Level   
        betaE  =  0.2;
   %--Curvature in the utility of Children's Expected Education  
        gamE    =  0.00;
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Education (Schooling) Levels
        EE_set= [0, 3.16, 6.98, 10.71, 14.68]';  % Using Brazilian data (1980)
                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Utility Cost of Living in the Different Regions
        tau_R = 0.0;
        tau_F = 0.273; 
        tau_C = 0.0;
        tau_P = 0.05;    
        tau_1_P = tau_P  ; %housing cost in terms of lost labor hours
  
        tau_1_F = 0; %utility cost in terms of curvatures only for favelas
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Production Parameters
        XA = 9.47;
        XM = 5.5; 
        XS = 1.990; 
    
        varphi     =  0.54;
        factor_wn  =  varphi^( varphi/(1-varphi) )*(1-varphi)*(XM/XS^varphi)^(1/(1-varphi));

%% Housing Prices: location l,   
%% Housing Cost Options:
%--- HousingCostOption=1 :  xiH_l    =   xiH0_l*(1+pop_l)^xiH1_l;
%--- HousingCostOption=2 :  xiH_l    =   xiH0_l*(pop_l)^xiH1_l;
%--- HousingCostOption=3 :  xiH_l    =   xiH0_l* exp( xiH1_l*(pop_l));

        HousingCostOption=1;

        xiH0_R=0;
        xiH0_F=0.34;
        xiH0_C=2.20;
    
        xiH1_R=0;
        xiH1_F=0.23;
        xiH1_C=0.3;
           

%% Initial State: cross-section education for the country
        pi_E = [0.2512 0.4550 0.1294 0.0955 0.0689]; %1980 25-40
   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      Computational Aspects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %% Generate or use previously generated labor market skills distributions
        GenerateSkills = 0;  %% 0: use previously generated; % 1: generate using Brazilian data

    %% Grid Sizes for Labor Market Skills
       Nxr = 200/2;
       Nxn = 300/2;

    %% Some (initial or guessed) distribution of the population
        Pop_data = [0.2874; 0.10; 0.3]; %1980 25-40
        Pop_R = Pop_data(1);
        Pop_F = Pop_data(2);
        Pop_P = Pop_data(3);
        Pop_C = 1 - Pop_R - Pop_F-Pop_P;
        PoP   = 1;

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      LABOR MARKET SKILLS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% --Previously generated to match cross-section distribution of earnings across occupations and education groups in Brazil
        NE=5;
                
        Moments_Skills_e = [   -1.2224   -7.3039    0.6669    2.7420   -0.0180
                                0.2326   -4.6167    0.7264    2.4773    0.1971
                                0.5676   -2.6523    0.7762    1.9710    0.3421
                                0.6784   -2.2963    0.9511    2.2052   -0.6787
                                1.2115    1.1570    0.8432    1.1728    0.4740];                  
        mu_r_e    = Moments_Skills_e(:,1);
        mu_n_e    = Moments_Skills_e(:,2);
        sigma_r_e = Moments_Skills_e(:,3);
        sigma_n_e = Moments_Skills_e(:,4);
        rho_rn_e  = Moments_Skills_e(:,5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      ESTIMATED PARAMETERS FOR THE INTERGENERATIONAL MOBILITY IN SKILL FORMATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    Call the function TransProb.m
%%    Now, with the option of using either U of C,F separated

%% Constructing the Implied Transition Functions
        paramsC =  [    0.3752    0.0602   -0.2782    0.0603   -1.4588    0.1549    2.2017   -0.8731
                        0.3452    0.2609   -0.5171    0.4215    0.6049    0.2326   -0.6450    0.3876
                        0.5579    0.0578    2.6716   -0.2078    3.5317   -0.1805    1.8001    0.0899
                        -0.0527    0.2080   -0.3673    0.3916   -0.7364    0.7135   -0.9283    0.7899
                        -2.9930    1.3062   -1.5812    1.2301   -1.2127    1.3262   -1.7434    1.5541];            

        paramsF =  [    -1.3926    0.3519   -4.3194    0.8462    1.7253   -2.2400    5.6594   -6.0559
                        0.3826    0.1879   -1.3421    0.4963   -0.8398    0.1813    0.7335   -0.9144
                        0.1229    0.1407   -0.2644    0.3590    0.1099    0.1168   -0.3958   -0.2198
                        0.3636    0.0828    0.8828    0.0198   -0.1174    0.4303   -0.3179    0.1135
                        -0.5405    1.0020   -0.7428    1.0815   -0.8126    1.2632   -0.5466    1.0847];

        paramsP =  [    -1.3926    0.3519   -4.3194    0.8462    1.7253   -2.2400    5.6594   -6.0559
                        0.3826    0.1879   -1.3421    0.4963   -0.8398    0.1813    0.7335   -0.9144
                        0.1229    0.1407   -0.2644    0.3590    0.1099    0.1168   -0.3958   -0.2198
                        0.3636    0.0828    0.8828    0.0198   -0.1174    0.4303   -0.3179    0.1135
                        -0.5405    1.0020   -0.7428    1.0815   -0.8126    1.2632   -0.5466    1.0847];
               
        paramsU =  [   -1.3447    0.3221   -4.0053    0.6585   -6.5037    0.9764   -0.2165   -0.5487
                        1.5671   -0.0174    0.7875    0.1079    1.9935   -0.1216    1.2530   -0.0760
                        2.1359   -0.0939    1.4595    0.0846    3.2487   -0.0719    1.8383    0.1454
                        1.0022    0.0308    2.5224   -0.0673    3.3225   -0.0551    2.8724    0.0567
                        2.9370   -0.0408    4.1377   -0.0805    4.3052    0.0418    3.8181    0.2640];

        paramsR =   [   -1.4031    0.5588   -3.7842    0.7317 -153.0621    0.4551 -343.9319    0.6787
                        0.2595    0.2825   -2.5330    1.0551    3.1706   -3.5829  -11.6147    2.9642
                        -5.0708    1.8372   -2.3179    0.9077    0.1668   -0.1033   -1.9433    0.5911
                        2.4999   -0.0915    1.4746    0.1801    3.4392   -0.4036    1.8471    0.2782
                        0.0341   -0.0626    2.8190   -0.9117    3.5031   -1.0195   -0.9551    0.6864];


%% Grid for numerical integration
        uB_xr        =    8; % 2*max(z_b*sigma_r_e + mu_r_e);
        uB_xn        =    8; % 2*max(z_b*sigma_n_e + mu_n_e);
        lB_xr        =   -10; %2*min(-z_b*sigma_r_e + mu_r_e);
        lB_xn        =   -10; %2*min(-z_b*sigma_r_e + mu_r_e);

        xr=linspace(lB_xr,uB_xr,Nxr);
        xn=linspace(lB_xn,uB_xn,Nxn)';
        ones_xr=ones(size(xr));
        ones_xn=ones(size(xn));

        xrxr=kron(xr,ones_xn);
        xnxn=kron(xn,ones_xr);

        Xr=exp(xr);
        XrXr=exp(xrxr);
        Xn=exp(xn);
        XnXn=exp(xnxn);


%% A Bivariate Log-normal Skills Distribution for each e \in E 
    for e=1:NE,    
        mu_xr    =   mu_r_e(e);
        mu_xn    =   mu_n_e(e);
        sigma_xr =   sigma_r_e(e);
        sigma_xn =   sigma_n_e(e);
        Sigma_xr =   sigma_xr^2;
        Sigma_xn =   sigma_xn^2;
        rho_rn   =   rho_rn_e(e);
        cov_xrxn =   rho_rn*sigma_xr*sigma_xn;

        frn_e(:,:,e)       = (1/(2*pi*(1-rho_rn^2)^.5*sigma_xr*sigma_xn))*exp( -1/(2*(1-rho_rn^2))*( (xrxr-mu_xr).^2/Sigma_xr + (xnxn-mu_xn).^2/Sigma_xn   -2*rho_rn/sigma_xr/sigma_xn*(xrxr-mu_xr).*(xnxn-mu_xn) ));
        fXrXn_e(:,:,e)     = frn_e(:,:,e)./XrXr./XnXn;  % densities for the exponentiated variables
        fXr_e(:,:,e)   = frn_e(:,:,e)./XnXn;
        fXn_e(:,:,e)   = frn_e(:,:,e)./XrXr;

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parameters of the Big Loop (over locations)    
   dist    =   10;
   J       =   1;
   Jmax    =   100;
   rlx     =   0.2;
   
%% New Parameters for the Extension
   xiH1_P = 0.34;  
   tau_P  = 0;
   xiH0_P = 0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Big Loop: Computing the Equilibrium Locations and Equilibrium Prices, Quantities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%% Location Choices: An Initial guessed value to start the iterations
    q_R0  = Pop_R*ones(Nxn,Nxr,NE);
    q_F0  = Pop_F*ones(Nxn,Nxr,NE);
    q_P0  = Pop_F*ones(Nxn,Nxr,NE);
    q_C0  = Pop_C*ones(Nxn,Nxr,NE);
    q_U0  = q_F0 + q_C0 + q_P0;

    q_R  = q_R0;
    q_F  = q_F0;
    q_P  = q_P0;
    q_C  = q_C0;
    q_U  = q_F0 + q_C0;

    
%% ---Speeding Up the Iterations, as the loop closes in the convergence
%  ---use as the initial values in the fminsearch functions.
    p_SR0 = alphaS;
    p_SU0 = alphaS/alphaA;    
    p_M0  = alphaM/alphaA;    

    p_SR  = p_SR0;
    p_SU  = p_SU0;
    p_M  = p_M0;
    
while J<Jmax;    
   
%% Comment: I added the term 1/PoP to make sure that the population adds to one; numerical integrals only get it close to 1. I only cor it once, at the first iteration.     
%% Constructing the implied skill distribution across regions, given the country's education distribution, piE, the location choices, q_l, and the conditional skill distribution, frn_e(:,:,e) 
if NE==5,
    frn_R   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_R(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_R(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_R(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_R(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_R(:,:,5) );
    frn_F   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_F(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_F(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_F(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_F(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_F(:,:,5) );
    frn_P   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_P(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_P(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_P(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_P(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_P(:,:,5) );
    frn_C   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_C(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_C(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_C(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_C(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_C(:,:,5) );
elseif NE==4,
    frn_R   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_R(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_R(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_R(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_R(:,:,4));
    frn_F   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_F(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_F(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_F(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_F(:,:,4));
    frn_P   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_P(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_P(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_P(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_P(:,:,4));    
    frn_C   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_C(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_C(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_C(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_C(:,:,4));
end

%% Constructing the implied skill distribution across regions, given the country's education distribution, piE, the location choices, q_l, and the conditional skill distribution, frn_e(:,:,e) 
for ie=1:NE,
        Pi_E_R(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_R(:,:,ie) ) );
        Pi_E_F(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_F(:,:,ie) ) );
        Pi_E_P(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_P(:,:,ie) ) );        
        Pi_E_C(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_C(:,:,ie) ) );
end

% Normalizaing to one, i.e. cross-section education distribution in each region
    pi_E_R =Pi_E_R./sum(Pi_E_R);
    pi_E_F =Pi_E_F./sum(Pi_E_F);
    pi_E_P =Pi_E_P./sum(Pi_E_P);
    pi_E_C =Pi_E_C./sum(Pi_E_C);
    
% Average Education in each Region
    EE_R = pi_E_R*EE_set;
    EE_F = pi_E_F*EE_set;
    EE_P = pi_E_P*EE_set;
    EE_C = pi_E_C*EE_set;
    
% Construction the Conditional Probabilities      
    for ie0=1:NE,
            % Conditional Transition Probabilities
                Pr_ee_R(ie0,:) = TransProb(ie0,EE_R,paramsR)';
                Pr_ee_F(ie0,:) = TransProb(ie0,EE_F,paramsF)';
                Pr_ee_P(ie0,:) = TransProb(ie0,EE_F,paramsP)';
                Pr_ee_C(ie0,:) = TransProb(ie0,EE_C,paramsC)';
            % Conditional Expectation of Children's Education, for each household, by Region
                CEE_R(ie0) = Pr_ee_R(ie0,:)*EE_set;
                CEE_F(ie0) = Pr_ee_F(ie0,:)*EE_set;
                CEE_P(ie0) = Pr_ee_P(ie0,:)*EE_set;
                CEE_C(ie0) = Pr_ee_C(ie0,:)*EE_set;
    end
    

%% Some other auxiliary distributions   
    f_XrXn_R      = frn_R./XrXr./XnXn;
    Xr_f_XrXn_R   = frn_R./XnXn;
    Xn_f_XrXn_R   = frn_R./XrXr;

    f_XrXn_F      = frn_F./XrXr./XnXn;
    Xr_f_XrXn_F   = frn_F./XnXn;
    Xn_f_XrXn_F   = frn_F./XrXr;
    
    f_XrXn_P      = frn_P./XrXr./XnXn;
    Xr_f_XrXn_P   = frn_P./XnXn;
    Xn_f_XrXn_P   = frn_P./XrXr;

    f_XrXn_C      = frn_C./XrXr./XnXn;
    Xr_f_XrXn_C   = frn_C./XnXn;
    Xn_f_XrXn_C   = frn_C./XrXr;
    
    f_XrXn_U      =    f_XrXn_C +    f_XrXn_F+    f_XrXn_P;
    Xr_f_XrXn_U   = Xr_f_XrXn_C + Xr_f_XrXn_F+ Xr_f_XrXn_P;
    Xn_f_XrXn_U   = Xn_f_XrXn_C + Xn_f_XrXn_F+ Xn_f_XrXn_P;
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Equilibrium, conditional on Location Decisions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    wu=XA; % pA=1, numeraire


%% Rural Regions
%% Setting up the function that prices services and allocates labor in Rural Areas
% Populations
    pop_R                   =    trapz( Xr,trapz(Xn,   f_XrXn_R ));  % employment in routine tasks, rural areas

% calling the function to find the equilibrium in Rural Regions
    p_SR0=p_SR; % After first using an arbitrary initial value, just use the previous equilibrium values 

    fun0=@(rural_p) RuralEquilibriumFunction(rural_p,pop_R, XS,XA,alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_R);    
            options = optimset('TolFun',1e-8);
            options = optimset('MaxIter',2000);
            options = optimset('MaxFunEvals',5000);
    [p_SR,fval_rural,exitflag_Rural,output] = fminsearch(fun0,p_SR0,options);

%% Just Making Sure it is positive:
    p_SR=max(p_SR,0.001);

%% Equilibrium in Rural Areas: Implied wages and labor market earnings
    wrR                     =    p_SR*XS;
    ybar_R                  =    p_SR*cSbar*(1-alphaS)/alphaS+cAbar;       % income threshold for the consumption of Services (presumes zero housing costs in rural areas)
% Labor Markets
    yR                      =    max(wu,wrR*XrXr);
    IOCC_r_R                =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_R(wu<wrR*XrXr)   =    1;
    ICS_R                   =    zeros(size(XrXr));  % Indicator, positive 
    ICS_R(yR>=ybar_R)       =    1;
    LrR                     =    trapz( Xr,trapz(Xn, Xr_f_XrXn_R .*IOCC_r_R));  % effective labor in routine tasks, rural areas
    ErR                     =    trapz( Xr,trapz(Xn,    f_XrXn_R .*IOCC_r_R));  % employment in routine tasks, rural areas
    QA                      =    XA*(pop_R-ErR);
    QS_R                    =    XS*LrR;
% Consumption of Services
    yTop_R                  =    trapz( Xr,trapz(Xn, yR.*f_XrXn_R.*ICS_R));
    fTop_R                  =    trapz( Xr,trapz(Xn,     f_XrXn_R.*ICS_R));
% Equilibrium Pricing Condition
    p_SR1                  =    alphaS*(yTop_R-cAbar*fTop_R)/(QS_R+(1-alphaS)*cSbar*fTop_R);
    QA_eq           =   QA;
    PSR_eq          =   p_SR;


%% Urban Regions   
%% Finding the No-Services equilibrium conditions in Urban regions
    rw_rn0 =  varphi/(1-varphi);
    funM=@(wr_wn) M_OnlyFunction(wr_wn,varphi,Xr,Xn,XrXr,XnXn, f_XrXn_U, Xr_f_XrXn_U, Xn_f_XrXn_U );
        options = optimset('TolFun',0.00000000000001);
        options = optimset('MaxIter',1000);
        options = optimset('MaxFunEvals',4000);
    [rw_rn,fval_Monly,exitflag_M_Only,output] = fminsearch(funM,rw_rn0,options);

%% just making sure it is positive....
    rw_rn=max(0.001,rw_rn);
    
%% M only allocations
    IOCC_r_U_M_only                        =  zeros(size(XrXr));
    IOCC_n_U_M_only                        =  zeros(size(XrXr));
    IOCC_r_U_M_only(rw_rn*XrXr >  XnXn)    =  1;
    IOCC_n_U_M_only(rw_rn*XrXr <  XnXn)    =  1;
    L_r_M_only                      =  trapz( Xr,trapz(Xn,  Xr_f_XrXn_U.*IOCC_r_U_M_only));
    Emp_r_U_M_only                  =  trapz( Xr,trapz(Xn,     f_XrXn_U.*IOCC_r_U_M_only));
    L_n_M_only                      =  trapz( Xr,trapz(Xn,  Xn_f_XrXn_U.*IOCC_n_U_M_only));
    Emp_n_U_M_only                  =  trapz( Xr,trapz(Xn,     f_XrXn_U.*IOCC_n_U_M_only));
    wrn_M_only                      =  rw_rn;
    Q_M_only                        =  XM*(L_r_M_only)^varphi*(L_n_M_only)^(1-varphi);
    pSM_M_only                      =  (XM/XS)*varphi*(L_n_M_only/L_r_M_only)^(1-varphi);

%% Setting up the function that solves for Urban Equilibrium (prices and allocations)
% Populations
    pop_C                   =    trapz( Xr,trapz(Xn,   f_XrXn_C ));  % employment in routine tasks, rural areas
    pop_F                   =    trapz( Xr,trapz(Xn,   f_XrXn_F ));  % employment in routine tasks, rural areas
    pop_P                   =    trapz( Xr,trapz(Xn,   f_XrXn_P ));  % employment in routine tasks, rural areas    
    pop_U                   =    pop_C + pop_F;
    pop = [pop_R; pop_F; pop_P; pop_C];
    pop = pop./sum(pop);
    pi_E_U = (pi_E_F.*pop_F+pi_E_P.*pop_P+pi_E_C.*pop_C)./pop_U;

    %% Hosing Costs: in units of M goods required to construct a house
      xiH_R    =   0; 
    if HousingCostOption==1,
        xiH_F    =   xiH0_F*(1+pop_F)^xiH1_F;
        xiH_P    =   xiH0_P*(1+pop_P)^xiH1_P;
        xiH_C    =   xiH0_C*(1+pop_C)^xiH1_C;
    elseif HousingCostOption==2,
        xiH_F    =   xiH0_F*(pop_F)^xiH1_F;
        xiH_P    =   xiH0_P*(pop_P)^xiH1_P;
        xiH_C    =   xiH0_C*(pop_C)^xiH1_C;
    elseif HousingCostOption==3,
        xiH_F    =   xiH0_F*exp(xiH1_F * pop_F);
        xiH_P    =   xiH0_P*exp(xiH1_P * pop_P);
        xiH_C    =   xiH0_C*exp(xiH1_C * pop_C);
                
    end
     
%% Calling the Function for Urban Equilibrium 
% Using equilibrium prices from previous location's iteration as initial condition 
    p_SU0 = p_SU;    
    p_M0  = p_M;
    pp0=[p_SU0,p_M0]; 
% calling the function and setting its options 
    fun1=@(urb_p) UrbanEquilibriumFunction_poorneigh(urb_p,pop_C,pop_F,pop_P,xiH_F,xiH_P,xiH_C,XS,XM,varphi,QA_eq, alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_C,f_XrXn_F,f_XrXn_P, L_r_M_only,L_n_M_only,Q_M_only,pSM_M_only);
            options = optimset('TolFun',1e-9);
            options = optimset('MaxIter',2000);
            options = optimset('MaxFunEvals',5000);
            options = optimset('TolX',1e-10);
    [urban_prices,fval_urb1,exitflag_Urban,output] = fminsearch(fun1,pp0,options);

%% Recovering the Equilibrium    
    p_SU    =  max(0.001,  urban_prices(1));
    p_M     =   max(0.001, urban_prices(2));
  
%% constructing measures of skills distributions in Urban Regions (U= C + F)    
   f_XrXn_U   =  f_XrXn_C +f_XrXn_F+f_XrXn_P;  
    Xn_f_XrXn_C   =  XnXn.*f_XrXn_C;
    Xn_f_XrXn_F   =  XnXn.*f_XrXn_F;
    Xn_f_XrXn_P   =  XnXn.*f_XrXn_P;
    Xn_f_XrXn_U   =  XnXn.*f_XrXn_U;
    Xr_f_XrXn_C   =  XrXr.*f_XrXn_C;
    Xr_f_XrXn_F   =  XrXr.*f_XrXn_F;
    Xr_f_XrXn_P   =  XrXr.*f_XrXn_P;
    Xr_f_XrXn_U   =  XrXr.*f_XrXn_U;

    if p_SU/p_M < pSM_M_only,
        wrU     =       varphi*p_M*XM*(L_n_M_only/L_r_M_only)^(1-varphi);
        wnU     =   (1-varphi)*p_M*XM*(L_n_M_only/L_r_M_only)^(-varphi);
    else
        wrU     =    p_SU*XS;
        wnU     =    factor_wn*( p_M /(p_SU)^varphi )^(1/(1-varphi));  % these are the wages implied by the case when both S and M are produced at U
    end

%% Individual Earnings
    yU                              =    max(wrU*XrXr, wnU*XnXn);   

%% Consumption of Services
    ybar_F                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_F;        % income threshold for the consumption of Services, households in slums
    ybar_P                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_P;        % income threshold for the consumption of Services, households in slums
    ybar_C                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_C;        % income threshold for the consumption of Services, households in cities
    ICS_F                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in slums
    ICS_P                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in slums
    ICS_C                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in cities
    ICS_F(yU>=ybar_F)               =    1;
    ICS_P(yU>=ybar_F)               =    1;
    ICS_C(yU>=ybar_C)               =    1;
% Aggregation: Consumption of Services
    yTop_F                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_F.*ICS_F));    % Slums
    fTop_F                  =    trapz( Xr,trapz(Xn,     f_XrXn_F.*ICS_F));
    yTop_P                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_P.*ICS_P));    % Poor neighborhood
    fTop_P                  =    trapz( Xr,trapz(Xn,     f_XrXn_P.*ICS_P));
    yTop_C                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_C.*ICS_C));    % Cities
    fTop_C                  =    trapz( Xr,trapz(Xn,     f_XrXn_C.*ICS_C));

% Aggregation: Supplies of Skills and Employment    
    IOCC_r_U                        =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_U(wnU*XnXn<wrU*XrXr)     =    1;
    LrU                             =    trapz( Xr,trapz(Xn, Xr_f_XrXn_U .*IOCC_r_U));      %  effective labor in routine tasks, urban areas
    ErU                             =    trapz( Xr,trapz(Xn,    f_XrXn_U .*IOCC_r_U));      %  employment in routine tasks, urban areas
    LnU                             =    trapz( Xr,trapz(Xn, Xn_f_XrXn_U .*(1-IOCC_r_U)));  %  effective labor in non-routine tasks, urban areas
    EnU                             =    trapz( Xr,trapz(Xn,    f_XrXn_U .*(1-IOCC_r_U)));  %  employment in non-routine tasks, urban areas
    if p_SU/p_M<pSM_M_only,
        QM                              =   Q_M_only;
        QSU                             =   0; 
    else
        LrM                                       =   (varphi*XM*p_M/XS/p_SU)^(1/(1-varphi))*LnU;
        QM  =   XM*(LrM)^varphi*(LnU)^(1-varphi);
        QSU =   XS*(LrU-LrM);          
    end

    p_SU1  =    alphaS*( yTop_F+yTop_P+yTop_C - fTop_C*(cAbar+p_M*xiH_C) - fTop_F*( cAbar+p_M*xiH_F ) - fTop_P*( cAbar+p_M*xiH_P ) )/( QSU +(1-alphaS)*cSbar*(fTop_F+fTop_P+fTop_C) );
    p_M1  =    alphaM/alphaA*( QA_eq - cAbar )/( QM -xiH_C*pop_C-xiH_F*pop_F-xiH_P*pop_P );
   
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%    Attained Utilities, Average Education Levels per Region, Location Choices
%% Consumption
    %% Indirect utilities: consumption for each region and skills levels
        Factor_Below  = (cSbar)^alphaS*(alphaA)^alphaA*(alphaM)^alphaM./(alphaA+alphaM)^(alphaA+alphaM);
        Factor_Above  = (alphaS)^alphaS*(alphaA)^alphaA*(alphaM)^alphaM;

    %% Indirect Utilities of Consumption: Baseline: low value assuming non-feasibility
        IndU_Cons_R=-9999*ones(Nxn,Nxr);
        IndU_Cons_F=-9999*ones(Nxn,Nxr);
        IndU_Cons_P=-9999*ones(Nxn,Nxr);
        IndU_Cons_C=-9999*ones(Nxn,Nxr);

    %% Minimum Income to Subsist in each Region
        y_subs_R = cAbar;
        y_subs_F = cAbar + xiH_F*p_M;
        y_subs_P = cAbar + xiH_P*p_M;
        y_subs_C = cAbar + xiH_C*p_M;

    if gam==1,
    %% Utility in Each Region, assuming no consumption of market services
        poor_IndU_Cons_R = log( (1-tau_R)*Factor_Below*(yR-y_subs_R)./(p_M)^alphaM );
        poor_IndU_Cons_F = log( (1-tau_F)*(Factor_Below^(1-tau_1_F))*((yU-y_subs_F).^(1-tau_1_F))./(p_M)^(alphaM*(1-tau_1_F)) );
        poor_IndU_Cons_P = log( (1-tau_P)*Factor_Below*((1-tau_1_P)*yU-y_subs_P)./(p_M)^alphaM );
        poor_IndU_Cons_C = log( (1-tau_C)*Factor_Below*(yU-y_subs_C)./(p_M)^alphaM );

    %% Utility in Each Region, assuming positive consumption of market services, {ICS_R, ICS_F, ICS_C=1}
        rich_IndU_Cons_R = log( (1-tau_R)*Factor_Above*(yR + p_SR*cSbar - y_subs_R)./( (p_M)^alphaM*(p_SR)^alphaS ) );
        rich_IndU_Cons_F = log( (1-tau_F)*(Factor_Above^(1-tau_1_F))*((yU + p_SU*cSbar - y_subs_F).^(1-tau_1_F))./( ((p_M)^alphaM*(p_SU)^alphaS )^(1-tau_1_F)) );
        rich_IndU_Cons_P = log( (1-tau_P)*Factor_Above*((1-tau_1_P)*yU + p_SU*cSbar - y_subs_P)./( (p_M)^alphaM*(p_SU)^alphaS ) );    
        rich_IndU_Cons_C = log( (1-tau_C)*Factor_Above*(yU + p_SU*cSbar - y_subs_C)./( (p_M)^alphaM*(p_SU)^alphaS ) );    
    else
    %% Utility in Each Region, assuming no consumption of market services
        poor_IndU_Cons_R = 1/(1-gam)*( (1-tau_R)*Factor_Below*(yR-y_subs_R)./(p_M)^alphaM ).^(1-gam);
        poor_IndU_Cons_F = 1/(1-gam)*( (1-tau_F)*(Factor_Below^(1-tau_1_F))*((yU-y_subs_F).^(1-tau_1_F))./((p_M)^alphaM^(1-tau_1_F)) ).^(1-gam);
        poor_IndU_Cons_P = 1/(1-gam)*( (1-tau_P)*Factor_Below*((1-tau_1_P)*yU-y_subs_C)./(p_M)^alphaM ).^(1-gam);
        poor_IndU_Cons_C = 1/(1-gam)*( (1-tau_C)*Factor_Below*(yU-y_subs_C)./(p_M)^alphaM ).^(1-gam);

    %% Utility in Each Region, assuming positive consumption of market services, {ICS_R, ICS_F, ICS_C=1}
        rich_IndU_Cons_R = 1/(1-gam)*( (1-tau_R)*Factor_Above*(yR + p_SR*cSbar - y_subs_R)./( (p_M)^alphaM*(p_SR)^alphaS ) ).^(1-gam);
        rich_IndU_Cons_F = 1/(1-gam)*( (1-tau_F)*(Factor_Above^(1-tau_1_F))*((yU + p_SU*cSbar - y_subs_F).^(1-tau_1_F))./( ((p_M)^alphaM*(p_SU)^alphaS )^(1-tau_1_F)) ).^(1-gam);
        rich_IndU_Cons_P = 1/(1-gam)*( (1-tau_C)*Factor_Above*((1-tau_1_P)*yU + p_SU*cSbar - y_subs_C)./( (p_M)^alphaM*(p_SU)^alphaS ) ).^(1-gam);
        rich_IndU_Cons_C = 1/(1-gam)*( (1-tau_C)*Factor_Above*(yU + p_SU*cSbar - y_subs_C)./( (p_M)^alphaM*(p_SU)^alphaS ) ).^(1-gam);
    end

    %% Attained Utilities: Consumption 
        %% First layer: Can the household subsist?   
            IndU_Cons_R(yR>y_subs_R) = poor_IndU_Cons_R(yR>y_subs_R);
            IndU_Cons_F(yU>y_subs_F) = poor_IndU_Cons_F(yU>y_subs_F);
            IndU_Cons_P(yU>y_subs_P) = poor_IndU_Cons_F(yU>y_subs_P);
            IndU_Cons_C(yU>y_subs_C) = poor_IndU_Cons_C(yU>y_subs_C);

        %% Second layer: Would the household buy Services in the market?
            IndU_Cons_R(ICS_R==1) = rich_IndU_Cons_R(ICS_R==1);
            IndU_Cons_F(ICS_F==1) = rich_IndU_Cons_F(ICS_F==1);
            IndU_Cons_P(ICS_F==1) = rich_IndU_Cons_P(ICS_F==1);
            IndU_Cons_C(ICS_C==1) = rich_IndU_Cons_C(ICS_C==1);

%% Children's Expected Education (CEE)
    %% Utility for CEE for each region and education levels     
    if gamE==1,
        for ie=1:NE,
            WCEE_R(ie) = log(CEE_R(ie));
            WCEE_F(ie) = log((1-tau_F)*(CEE_F(ie)^(1-tau_1_F)));
            WCEE_P(ie) = log(CEE_P(ie));
            WCEE_C(ie) = log(CEE_C(ie));
        end
    else
        for ie=1:NE,
            WCEE_R(ie) = 1/(1-gamE)*(CEE_R(ie))^(1-gamE);
            WCEE_F(ie) = 1/(1-gamE)*((1-tau_F)*(CEE_F(ie)^(1-tau_1_F)))^(1-gamE);
            WCEE_P(ie) = 1/(1-gamE)*(CEE_P(ie))^(1-gamE);
            WCEE_C(ie) = 1/(1-gamE)*(CEE_C(ie))^(1-gamE);
        end
    end
   
  %% Attained Utilities: U[Consumption] + betaE*CEE_l
  %    
  for ie=1:NE,
        V_R(:,:,ie) = (1-betaE)*IndU_Cons_R + betaE*WCEE_R(ie);
        V_F(:,:,ie) = (1-betaE)*IndU_Cons_F + betaE*WCEE_F(ie);
        V_P(:,:,ie) = (1-betaE)*IndU_Cons_P + betaE*WCEE_P(ie);
        V_C(:,:,ie) = (1-betaE)*IndU_Cons_C + betaE*WCEE_C(ie);
  end

%% Equilibrium Location Decisions
    q_R1 = exp(etaE_R*V_R)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_P) + exp(etaE*V_C));
    q_F1 = exp(etaE*V_F)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_P) + exp(etaE*V_C));
    q_P1 = exp(etaE*V_P)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_P) + exp(etaE*V_C));
    q_C1 = exp(etaE*V_C)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_P) + exp(etaE*V_C));

    dist          =  max(max(max(abs(q_C-q_C1))))+ max(max(max(abs(q_P-q_P1))))+ max(max(max(abs(q_F-q_F1))))+ max(max(max(abs(q_R-q_R1))));
    Dist_store(J) = dist;
    
    q_R= ( rlx*q_R1+(1-rlx)*q_R );
    q_F= ( rlx*q_F1+(1-rlx)*q_F );
    q_P= ( rlx*q_P1+(1-rlx)*q_P );
    q_C= ( rlx*q_C1+(1-rlx)*q_C );
    
    pi_E_L = [pi_E_R; pi_E_F; pi_E_P; pi_E_C]';
    pi_E_calc = pi_E_L*pop;
    pi_E_calc = pi_E_calc./sum(pi_E_calc);
    pi_E = ( rlx*pi_E_calc'+(1-rlx)*pi_E ); 
    
    disp([' Iteration: ',num2str(J),' distance: ',num2str(dist)])
    J=J+1;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reporting!

%% Total Output, 
    GDP    =  QA+p_M*QM+p_SR*QS_R+p_SU*QSU;
    VA_A   = 100*QA/GDP;
    VA_SR  = 100*p_SR*QS_R/GDP;
    VA_SU  = 100*p_SU*QSU/GDP;
    VA_M   = 100*p_M*QM/GDP;
    
    GNI    =  wu*(pop_R-ErR)+wrR*LrR+wrU*LrU+wnU*LnU;
    
    PYS_u  = 100*wu*(pop_R-ErR)/GNI; 
    PYS_rR  = 100*wrR*LrR/GNI; 
    PYS_rU  = 100*wrU*LrU/GNI; 
    PYS_nU  = 100*wnU*LnU/GNI; 

format short


%% Results

     l_A = (pop_R-ErR);
     l_SR = ErR;
     l_SU = (1-LrM/LrU)*ErU;
     l_S = l_SR + l_SU;
     l_M = (LrM/LrU)*ErU + EnU;
     
     EE_U = (EE_F*(pop_F/(pop_F+pop_P+pop_C))+EE_P*(pop_P/(pop_F+pop_P+pop_C))+EE_C*(pop_C/(pop_F+pop_P+pop_C)));
     
     EE_R_future = pi_E_R*CEE_R';
     EE_F_future = pi_E_F*CEE_F';
     EE_P_future = pi_E_P*CEE_P';
     EE_C_future = pi_E_C*CEE_C';
     EE_U_future = (EE_F_future*(pop_F/(pop_F+pop_P+pop_C))+EE_P_future*(pop_P/(pop_F+pop_P+pop_C))+EE_C_future*(pop_C/(pop_F+pop_P+pop_C)));
     
     pi_E_comp(1)=pi_E_R(1)*pop_R+pi_E_F(1)*pop_F+pi_E_P(1)*pop_P+pi_E_C(1)*pop_C;
     pi_E_comp(2)=pi_E_R(2)*pop_R+pi_E_F(2)*pop_F+pi_E_P(2)*pop_P+pi_E_C(2)*pop_C;
     pi_E_comp(3)=pi_E_R(3)*pop_R+pi_E_F(3)*pop_F+pi_E_P(3)*pop_P+pi_E_C(3)*pop_C;
     pi_E_comp(4)=pi_E_R(4)*pop_R+pi_E_F(4)*pop_F+pi_E_P(4)*pop_P+pi_E_C(4)*pop_C;
     pi_E_comp(5)=pi_E_R(5)*pop_R+pi_E_F(5)*pop_F+pi_E_P(5)*pop_P+pi_E_C(5)*pop_C;
     
     pi_E_R_future(1) = pi_E_R*Pr_ee_R(:,1);
     pi_E_R_future(2) = pi_E_R*Pr_ee_R(:,2);
     pi_E_R_future(3) = pi_E_R*Pr_ee_R(:,3);
     pi_E_R_future(4) = pi_E_R*Pr_ee_R(:,4);
     pi_E_R_future(5) = pi_E_R*Pr_ee_R(:,5);
     
     pi_E_F_future(1) = pi_E_F*Pr_ee_F(:,1);
     pi_E_F_future(2) = pi_E_F*Pr_ee_F(:,2);
     pi_E_F_future(3) = pi_E_F*Pr_ee_F(:,3);
     pi_E_F_future(4) = pi_E_F*Pr_ee_F(:,4);
     pi_E_F_future(5) = pi_E_F*Pr_ee_F(:,5);  
     
     pi_E_P_future(1) = pi_E_P*Pr_ee_P(:,1);
     pi_E_P_future(2) = pi_E_P*Pr_ee_P(:,2);
     pi_E_P_future(3) = pi_E_P*Pr_ee_P(:,3);
     pi_E_P_future(4) = pi_E_P*Pr_ee_P(:,4);
     pi_E_P_future(5) = pi_E_P*Pr_ee_P(:,5);     
     
     pi_E_C_future(1) = pi_E_C*Pr_ee_C(:,1);
     pi_E_C_future(2) = pi_E_C*Pr_ee_C(:,2);
     pi_E_C_future(3) = pi_E_C*Pr_ee_C(:,3);
     pi_E_C_future(4) = pi_E_C*Pr_ee_C(:,4);
     pi_E_C_future(5) = pi_E_C*Pr_ee_C(:,5); 
     
     pi_E_future(1) = pi_E_R_future(1)*pop_R+pi_E_F_future(1)*pop_F+pi_E_P_future(1)*pop_P+pi_E_C_future(1)*pop_C;
     pi_E_future(2) = pi_E_R_future(2)*pop_R+pi_E_F_future(2)*pop_F+pi_E_P_future(2)*pop_P+pi_E_C_future(2)*pop_C;
     pi_E_future(3) = pi_E_R_future(3)*pop_R+pi_E_F_future(3)*pop_F+pi_E_P_future(3)*pop_P+pi_E_C_future(3)*pop_C;
     pi_E_future(4) = pi_E_R_future(4)*pop_R+pi_E_F_future(4)*pop_F+pi_E_P_future(4)*pop_P+pi_E_C_future(4)*pop_C;
     pi_E_future(5) = pi_E_R_future(5)*pop_R+pi_E_F_future(5)*pop_F+pi_E_P_future(5)*pop_P+pi_E_C_future(5)*pop_C;
     
     C_M=(QM-xiH_F*pop_F-xiH_P*pop_P-xiH_C*pop_C);
     C=QA+p_M*C_M+p_SR*QS_R+p_SU*QSU;
     C_A_share=QA/C;
     C_S_share=(p_SR*QS_R+p_SU*QSU)/C;


results = [p_SR p_SU p_M wu wrR wrU wnU  0 xiH_F*p_M xiH_C*p_M ...
                QA_eq QS_R QSU QM GDP VA_A VA_SR VA_SU VA_SR+VA_SU VA_M GNI PYS_u PYS_rR PYS_rU PYS_nU...
                l_SR l_SU l_S l_M l_A pop_R pop_F pop_P pop_C...
                EE_R,EE_F,EE_P,EE_C EE_U ...
                EE_R_future,EE_F_future,EE_P_future,EE_C_future, EE_U_future ...
                pi_E_R(1) pi_E_R(2) pi_E_R(3) pi_E_R(4) pi_E_R(5) ...
                pi_E_F(1) pi_E_F(2) pi_E_F(3) pi_E_F(4) pi_E_F(5) ...
                pi_E_P(1) pi_E_P(2) pi_E_P(3) pi_E_P(4) pi_E_P(5) ...
                pi_E_C(1) pi_E_C(2) pi_E_C(3) pi_E_C(4) pi_E_C(5) ...
                pi_E_U(1) pi_E_U(2) pi_E_U(3) pi_E_U(4) pi_E_U(5) ...
                pi_E_future(1) pi_E_future(2) pi_E_future(3) pi_E_future(4) pi_E_future(5)];
            
 params = [xiH1_P tau_P xiH0_P tau_1_P  alphaS alphaA ...
           cSbar cAbar gamE etaE etaE_R betaE gamE ... 
           tau_R tau_F tau_C tau_1_F ...
           XA XM XS varphi ...
           xiH0_R xiH0_F xiH0_C xiH1_R xiH1_F xiH1_C];
            
    % fid = fopen( 'exceloutputvar.txt', 'a+' );
    % fprintf(fid, '%18.10f\n', results);
    % fclose(fid);       
    % 
    % fid = fopen( 'exceloutputparams.txt', 'a+' );
    % fprintf(fid, '%18.10f\n', params);
    % fclose(fid);

